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The ground state phase diagram of the extended Hubbard model containing nearest and next- 
to-nearest neighbor interactions is investigated in the thermodynamic limit using an exact method. 
It is found that taking into account local correlations and adding next-to-nearest neighbor inter- 
actions both have significant effects on the position of the phase boundaries. Improved stability 
domains for the ^-pairing state and for the fully saturated ferromagnetic state at half filling have 
been constructed. The results show that these states are the ground states for model Hamiltonians 
with realistic values of the interaction parameters. 

PACS number: 74.20.-z, 75.10.Jm, 75.10.Lp 

I. INTRODUCTION 

Exact solutions in physics are of great importance, since in some cases the errors introduced by the approximations 
may dominate the results to such an extent that one might end up with an incorrect description of the studied 
phenomenon. When applying an analytic, but non-exact, approach one has to know to what extent the approximation 
is valid. In case of perturbative methods a well-defined small parameter can ensure that the higher order terms are 
indeed negligible. Often, it is hard to find such a small parameter due to the fact that the given phenomenon itself 
has strong-coupling characteristics, i.e. the associated correlation effects are not small at all. This is especially true 
by the investigation of strongly correlated electron systems. Ferromagnetism is an example for intermediate-to- strong 
coupling phenomenon, where one has to be very cautious to apply perturbative approaches. 

Considering the exact results with respect to the dimensionality D, we can see that most of them have been derived 
in two limiting cases: either D = 1 or D = oo. For example, the exact solution of the Hubbard model was given in D = 1 
dimension by means of the Bethe-ansatz by Lieb and WuEl. The other class of exact solutions belongs to the other 
limiting case, i.e. D — oo, when the dynamical mean-field approximation becomes exactma. The situation, however, 
gets more complicated as physically interesting, lower dimensional cases (e.g. systems in D = 2 or D — 3 dimensions) 
are considered. D > 1 rules out the applicability of the well-established Bethe-ansatz approach, while mean-field like 
descriptions lead to qualitatively, or quantitatively incorrect conclusions, because the effects of spatial fluctuations are 
not taken properly into accountertl 

In recent years, some exact, non-perturbative methods have, been developed to investigate the ground state of 
Hubbard and Hubbard-like models in large parameter regimesou. In the presentj-work we focus on the so-called 
optimum ground state method (OGS) established by de Boer and SchadschneideiO. Using this method one can 
obtain rigorous constraints on the model parameters which define regions where e.g. the CDW, the Neel, the fully 
saturated ferromagnetic, or the ^-pairing state of momentum P becomes the exact ground state of the Hamiltonian. 

The basic idea of the OGS method is to diagonalize a specially chosen local Hamiltonian and to tune all the model- 
parameters such a way that all local eigenstates which are needed to the construction of a given global ground state 
are also local ground states. This means that, on one hand, the corresponding eigenvalues of the local Hamiltonians 
should be all equal in magnitude and, on the other hand, this common value should be the lowest eigenvalue of the 
local problem. Following this method, one can obtain different regions in the parameter space of the model defined 
by inequalities. These inequalities mean sufficient conditions for a state to be the ground state inside a special region. 
Outside the derived region the state under study may or may not be the ground state of the model. 

There are basically two different ways to enlarge the region of guaranteed stability: taking a larger local Hamiltonian 
or incorporating next-to-nearest neighbor interactions. As far as the first approach is concerned, it is obvious that 
the extent of local correlations which are taken into account is controlled by the size of the local Hamiltonian to be 
diagonalized exactly. Therefore, using local Hamiltonians defined on larger clusters of the lattice should tipically lead 
to better constraints (i.e. extended stability regions), even if purely nearest neighbor interactions are present. 

It is also well-knossai that next-to-nearest neighbor hopping has important effects. For instance, it was shown 
rigorously by Tasakiliil that the pure Hubbard model characterized by hopping of electrons between nearest and 
next-to-nearest neighboring sites with dispersive bands exhibits ferromagnetism for_£nite Coulomb interaction at 
zero temperature. Recent Projection Quantum Monte Carlo studies by Hlubina et al.tB confirmed this fact for finite 
temperatures, too. Beside the consequences of longer range hopping, the importance of nearest and next-to-nearest 
neighbor off-site interactions (both diagonal and off-diagonal) has also been emphasized both from experimental and 
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from theoreticalH113'tlHlZHl§ sides. These extra terms (density-density type interaction, correlated hopping of electrons, 
hopping of electron pairs and the exchange coupling) all originate from the spin independent Coulomb interaction of 
electrons. Nevertheless, the values of these longer range interactions (e.g. between next-to-nearest neighboring sites) 
are known neither theoretically nor experimentally because of the complicated nature of screening processes in solids. 
However, it is obvious that these interactions are present in real materials. Their strength decreases with increasing 
interatomic distances on the lattice and they can have important effects on the characteristics of strongly correlated 
electron systems. Hence, it is a challenging task to incorporate and to treat them in an exact way on the level of the 
model Hamiltonian. l_ . 

The aim of the present paper is two-fold. First, we would like to extend the previous calculations of Ref.EJ 
by choosing a larger local Hamiltonian which is defined on elementary plaquettes consisting of four lattice sites of 
the D dimensional hypercubic lattice. Using these local Hamiltonians and a simple, numerically exact method we 
have constructed the stability domains for the ry-pairing states of momentum P — 0, ir and for the fully saturated 
ferromagnetic state in the parameter space of an extended Hubbard model with a half-filled band. Given the size of 
the local Hamiltonian, the diagonalization is done numerically. The stability regions are deduced from the equality of 
the lowest eigenvalue of the chosen local Hamiltonian and of an upper bound derived appropriately from the variational 
principle of quantum mechanics. Second, the present choice of the larger local Hamiltonian gives us a simple way to 
incorporate next-to-nearest neighbor interactions, to treat them exactly and to investigate their effects on stability. 
The paper is organized as follows: in Sec. |l|we intr oduc e the meth od. In Sec. Ill the global Hamiltonian and the 



corresponding local Hamiltonian are defined. Sections IV A and [VB contain the stability domains in D — 2, 3 for the 



?7-pairing states of momentum P and for the fully saturated ferromagnetic case, respectively. Finally a short summary 
and discussion closes the presentation in Sec. M. 

II. METHOD 

Let us consider a model Hamiltonian defined on a discrete lattice. This (so-called global) Hamiltonian can be 
decomposed into the sum of equivalent local Hamiltonians h c i ustcl - defined on identical clusters of lattice sites, the 
union of which covers the lattice. In other words 

H = 2_j ^cluster • (1) 
all clusters 

If only on-site interactions are considered, each cluster might consist of a single lattice point and /i c iustcr is simply 
a Hamiltonian defined on the individual lattice points. If intersite interactions are also present, the cluster has to 
consist of (at least) two lattice sites (in case of only nearest neighbor interactions) and the proper minimal local 
Hamiltonian is actually a Hamiltonian defined on a bond joining two (in the simplest case nearest neighboring) sites. 
Nevertheless, the incorporation of longer and longer range interactions, or of more and more spatial correlations, 
requires the enlargement of the minimal cluster and hence the corresponding local Hamiltonian Cluster- In principle, 
even infinite range interactions and correlations can be taken into account, but this would require Cluster — H. The 
tractable cluster size is, however, limited by the feasibility of the necessary exact diagonalizations. 
Setting up the local Hamiltonian and choosing a suitable local basis, the eigenvalue problem 

^cluster | ^cluster) — ^cluster|0cluster) (2) 

can be solved exactly. Thus the full spectrum e*i uster (i = 1, • • ■ , dim(ft, c i ustcr )) of Cluster can be obtained. Since the 
clusters of decomposition are equivalent, the exact ground state energy Eqs is bounded from below by the relation 

-Elower = -^cluster mm Cluster — ^GS ! (3) 

i 

where the number of clusters on the considered lattice can be written as iV c i us t e r = fL. Here / represents a simple 
combinatorical factor depending on the dimension and structure of the underlying lattice and L stands for the number 
of lattice sites. 

The variational principle of quantum mechanics provides an upper bound for Eqs , namely 

„ , (*trial|£f|*trial) „ ,.s 

\* trial I Atrial/ 

where Atrial) stands for an arbitrary trial wavefunction. In our case Atrial) is an exact eigenstate of the global 
Hamiltonian for a certain set of model parameters. 
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Combining now Eqs. (§) and © yields 



Slower < SqS < Supper • (5) 

Exploiting that both Sower and S uppcr are analytic functions of the couplings of the global Hamiltonian, after carefully 
adjusting the coupling constants one can satisfy the equality 

Sl owcr = S U p pcr = SqS , (6) 

which means that for a certain set of model parameters (in a special sector of the ground state phase diagram) the 
exact ground state energy Eqs is found. Furthermore, if the ground state has no degeneracy, the exact ground state 
*gs) °f the global Hamiltonian is also found. In our case the non-degeneracy is provided by the fact that the states 



we consider (see Sec. IV) can be built up simply by using the lowest energy eigenstate of the local Hamiltonian. As an 
example let us consider the case of the D = 1 dimensional half-filled chain with bonds as clusters. For certain values of 
the coupling constants it can be reached that the non-degenerate lowest eigenvalue of the local Hamiltonian belongs 
to the parallel orientation of spins included on the bond. Since the bonds are equivalent, all the bonds along the chain 
contain parallel-oriented spins with the same local energy for the given set of model parameters. This fact yields the 
long range order of spins (here ferromagnetism) and also the non-degeneracy of the global ground state (apart from 
the spin-degeneracy). Similar arguments hold for the non-degeneracy of the ground state in higher dimensions, too. 

Changing now the trial wavefunction and following the procedure discussed above a new ground state in a different 
region of the phase diagram can be found. Furthermore, repeating the method with more and more trial wavefunctions 
a large portion of the ground state phase diagram of the model can be explored. 



III. GLOBAL AND LOCAL HAMILTONIAN 



Let us define now our global Hamiltonian onafl dimensional hypercubic lattice (D>1) in the following form 

Sgiob = UU + [— Ujtij + XijXij + VijVij + YijYij + JijJij — /i/t , (7) 



ho 



where 



i 

(7 

Xij — c \a C jv ( n i,-a + n j,-cj) 

(7 



Jij — — A^ y (S^Sj + S^S i ) + Afj S*Sj 



Here the fermion operators c\ a (cj CT ) create (annihilate) electrons with spin a in the single tight-binding Wannier 
orbital associated with site i. is the particle number operator of electrons with spin cr, and ni = ni\ + ny. 
Furthermore, the spin operators are given by = ct^Cj^, 1 S l i ~=ct,Cjf and S^ = ^(n^ — n^). The above Hamiltonian 
contains a term corresponding to the familiar Hubbard term of doubly occupied sites (U), the hopping of a single 
electron (characterized by Uj), a density dependent (or correlated) hopping term (Xij), an intersite density-density 
type interaction of electrons (14,-), a term describing the hopping of electron pairs (Yij), for Aj^ y = Afj = 1 a 
Heisenberg-type exchange interaction of spins (Jij) and a chemical potential term (/j,). We note that = and 

Afj = 1 represents an Ising-type coupling of electron spins, while the case of = 1 and A^ = corresponds to an 
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XY-type interaction. Jy > ( Jy < 0) means antiferromagnetie, (ferromagnetic) type of exchange. The relevance of 
the present model for real materials is discussed e.g. in Refs.u'E3. 

In the rest of the paper we investigate only couplings between on-site, nearest-neighbor and next-to-nearest neighbor 
electrons. Incorporating this restriction into Eq. (Q) and using the convention Ai—Aij (1 — 1,2 for i, j being nearest 
and ncxt-to-nearcst neighboring sites, respectively) for the intersite couplings, one can rewrite the global Hamiltonian 
as 



tfgiob = uJ2 - \ J (mi - ~J (8) 
2 f i i 

+ Y Y \ 2 Y [Mm,-* + n jt - a ) - U] (c\ a c ja + c\ a c ia ) + - Vi(m - l)( nj - 1) 
+ \yMA&M +c )Ai<*M) + l^JiS+S- +S+Sr) + J«5?S?| 



IJ.fi - Eq 



Here X^ij), means a summation over nearest neighboring (^ = 1) and next-to-nearest neighboring (1 — 2) sites. The 
...(.<;. = /s.f Y Ji and — AfJi are also used. In the special case of Af Y — Af — 1, however, the 



J; = Jxy = Jz^ notation will be kept for the sake of clarity. Eq represents a numerical constant which shifts the zero 
point of the energy scale. We remark that during the reformulation of Eq. (fj]) to Eq. (||) the chemical potential is 
also shifted by some constant. 

Because we have only nearest and next-to-nearest neighbor interactions, the fully symmetric minimal clusters 
are the two-dimensional elementary plaquettes of the D — 2 dimensional square lattice, as depicted in Fig. fil All 
D>2 dimensional hypercubic lattices can be covered with the elementary plaquettes, however, in those cases the D 
dimensional hypercubes would be the fully symmetric minimal clusters. 

Following the method of Sec. [0] the global Hamiltonian can be rewritten in terms of such plaquettes (J2[i j i m ] 
means a summation over the plaquettes) as 

-Hglob = E Kjlm (9) 

[i,j,l,m\ 



with 



lijlm = £ Y (" Q T - ^ (jM - i 



^ Y Y ( C L<r C <3^ + 4,<t Cq ^) ( n «.-<7 + n 0,-a) 



4A 



J£ , V 



h 

jr Y Y ( c <U c /^ + c l,* Ca ^) - j Y Y ( c i,aC/3,<x + 4> <a c a , c 

J1 ( a ,f3)eA N <y 32 (a,f3)eA NN <y 

^ ]T (n a -l)(n -l) + ^ J2 

1 1 (a,p)eA N 12 (a,f3)£A NN 

47" Y ( c U c U c pi c pi + c U c li Ca ^) + ~r Y { c U c li c 0i c 0i + c l\ c li Ca t Ca ' 

3l (a,f3)eA N J2 (a,f3)£A NN 

t(1) . , t(2) 



(a,/3)eA N (u,f3)£A NN 

+ w Y s « s l + J i 2) Y s « s p -jY n * ( 10 ) 

Jl (a,j3)eA N (a,f3)£A N N 2 a&Ao 

Here /i =22/^1 and J2 — 2 are numerical constants (z\ — 2D and Z2 = 4(^) are the number of nearest and next- to- 
nearest neighboring sites, respectively, on the D dimensional hypercubic lattice). They are needed to avoid double 
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or higher counting of intersite interactions during the plaquette summation. Furthermore, t\, X\, V\, Y\ are the 
values of single electron hopping, correlated hopping, density-density type interaction and pair-hopping between 
nearest neighboring sites, respectively, while t 2 , A 2 , V%, Y 2 indicate the analogous processes between next-to- nearest 
neighboring sites. U is the Hubbard interaction which can either be positive or negative in our model. To mimic 

real systems, however, it should be repulsive. In addition, spin interactions with exchange couplings jiy\ jf$ , 

(2) 

and Jz are included on the plaquette. Further notations: Ao = {i, j,l,m} means the set of individual lattice points, 
•AN = {(i,j), (ji 0: {h m )i ( m i i)} represents the set of nearest neighboring sites and ANN = {(h 0> (i> m )l indicates the 
set of next-to- nearest neighboring sites on the plaquette (see Fig. |l|). The possibility of anisotropies can be naturally 
incorporated into the local (and also global) Hamiltonian via the non-equivalence of the orthogonal directions of the 
plaquette. The effects, however, are not discussed in the present paper. 

To apply Eq. (^) the connection between the number of lattice sites L and the number of clusters iV c i UB t er (here 
plaquettes) is also needed. Combinatoric considerations give the following simple result for this quantity on a D 
dimensional hypercubic lattice, 

^cluster = \D{D - l)L . (11) 



IV. RESULTS 



As an illustration of the method described in Sec. |n| we consider a few physically interesting states, the ry-pairing 
states of momentum P (Sec. IV A) which s how o ff-diagonal long range order and are hence superconducting, and 
the fully saturated ferromagnetic state (Sec. IV B). We determine under what conditions these states are the ground 
states of the global Hamiltonian i? g i b- All the calculations presented are done at half filling. Except the repairing 
state of momentum P = 0, it is not possible to express the results in a compact analytic form as it has been done 
earlier in Ref.E3. Therefore, and for the sake of visualization, we have restricted ourselves to special cuts of the 
parameter space in illustrating the effects of the larger local Hamiltonians and of next-to-nearest neighbor couplings. 
The cuts are chosen in such a way that the corresponding ground state ]Afse diagrams can easihz-.be compared with 
the previously published rigorous results of Strack et al&H, de Boer et al£3'E3 and Montorsi et alEIl. In principle, one 
can also investigate the role of each nearest and next-to-nearest neighbor interactions separately, one by one, using 
our method. 



A. ^-pairing states of momentum P 



The definition of the ^-pairing operator of momentum P is given by the relation 

L 



4 = E 



Using this operator an repairing state of momentum P and of pairs N can be constructed as 

JV 



= K (4) |0), 



(12) 



(13) 



where K - 



(L-N)\ 

li m 



is a normalization factor. For further details about the ^-pairing states the reader is referred 

to the literature£^l~E3. Since we would like the ^-pairing states to be the ground state of our model, it is instructive 
to calculate the commutator of the ^-operator with the global Hamiltonian if g i b- One finds that 



- (X k t k ) J2 (e iPj + e lPl ) (4i c h + c U c 5t) 
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(ji)k k v 



t ,t 
u 



Ffee 



- t^*") ( nj - 1) 4 tC ; 

- U k e lP ^ (ni - 1) c}^ 



(14) 



Calculating now the quantity 



-ffgiob, (f/p) 



|0), one can easily deduce the parameters where the ?7-pairing states of 

momentum P are the ground state of the starting model; for momentum P = one arrives at the requirements Xi = ti 
and Yi=2Vi (£ = 1,2), while for momentum P = ir the conditions X 2 = t2 and Yj = (— l) s 2Vi (£ = 1,2) must be satisfied. 
One would also look for 77-pairing states of momentum P ^ or P =/= n. These, however, represent the ground states 
of tbG,rgJobal Hamiltonian in Eq. (0), where Xi=ti (£ = 1,2), U < —At\ and all the other interaction constants are 
zeroc3e3. 

Using the 77-pairing states as trial wavefunctions, the upper bound in the thermodynamic limit for the ground state 
energy per lattice site L is 



£"' p 

upper 



1 



{U + ziVi +Z2V2) + 



1 i l 
— n — n 

2 y 2 



2 

1=1 



ziVi + -n(l - ^n) ZjYi cos 1 P - fin 
1=1 



(15) 



Exploiting the constraints between the amplitudes of pair-hopping Yi and that of density-density type interaction Vi 
one gets for the half filled case (n=l) that 



=-L 

upper ^ 



u 



2 

E 

1=1 



ziVi - 11L 



(16) 



Despite the fact, that the upper bounds for the 77-pairing states of momentum P — and P — tt are identical, there is 
a characteristic difference between the two sets of wavefunctions; the state |^(N, P = it)) remains an exact eigenstate 
of the global Hamiltonian 7? g i b even if X\ ^ t\ . 

As mentioned earlier, it is possible to include more spatial correlations using local Hamiltonians defined on larger 
clusters of the lattice. This leads to the extension of the stability of the chosen state and means an improvement 
on the phase boundaries. For the ^-pairing state of momentum P = this effect is depicted in Fig. ||, where the 
J xy -J z cut of the coupling constants' space is chosen. The-inner triangle corresponds to the stability domain of the 
fyo-state determined by the OGS method of de Boer et alHj using bond Hamiltonians with purely nearest neighbor 
interactions. The shaded regions give the improvements of the boundaries applying our method with plaqucttc 
Hamiltonians containing only nearest neighbor interactions. The axes J a = J& A(U, V%) (a = xy,z) represent the 
rescaled values of nearest neighbor exchange interactions with the scaling factor of 



A(U,Vx) = 





2U J T 




2 


— + Vi 






Zl 





(17) 



From Fig. || one can immediately read off the stability criteria for the 770-state to be the ground state of -ffgiob ■ In the 
absence of next-to-nearest neighbor interactions (i.e. i2 = ^2 = ^2 = 12 = ^2 = 0) for Vi < values of nearest neighbor 
density-density type interaction 



V\ < 
-1 < ^ < 

~ h ~ 



J V 



(18) 







This is a considerable improvement over the V\ < 0, — 1 < J z /ti < —2 | J xy /ti + 1 criteria derived following Ref. 
using bond Hamiltonians. 

If now next-to- nearest neighbor interactions are turn on, i.e. nonzero values of V2 (and hence Yj) are considered, 
the criteria of Eq. ( |l8| ) remains valid in an improved form. In this case there is a need for the proper change of the 
scaling factor A, which now depends also explicitly on next-to-nearest neighbor interactions. The new value of the 
scaling factor is 



A(U,V U V 2 ) 



£2 

zi 



2U 

Z2 



+ Vi 



v 2 



(19) 



G 



In real systems the X\ = t\ requirement does not hold in general. However, the strengths of the corresponding 
two interactions are of the same magnitude. For X\ ^ t\ only the //-pairing state of momentum P = n represents 
an exact eigenstate of the model Hamiltonian i? g i b- Figure ^ shows the stability regions of the 7^-state for two 
different sets of model parameters expressed in units of eVin the U-t\ plane. Dotted lines represent the boundary for 
stability regions corresponding to bond Hamiltonians, while solid lines are the boundaries calculated with plaqucttc 
Hamiltonians in the absence of next-to-nearest neighbor interactions. The shaded regions clearly show the extension 
of stability regimes due to the choice of larger local Hamiltonians. 

We now consider the phase diagrams in the U-Y\ (Fig. |j) and V2-V1 (Fig. ||) planes for fixed values of the other 
interactions. In Fig. [I] the exchange interaction has been fixed by the relation Ji = — 2Y±. This assures that the 



global Hamiltonian of Sec. |IH| copcides with the model Hamiltonian of Ref.tfj, where the authors, using the method 1 
positive semi-definite operatorsH'E30, derived rigorous bounds for the Ti^-state. Comparing Fig. ^ with Fig. [l]in Ref.tJ 
two basic differences can be noticed. First, the boundary of the stability region, even for the special case of Xx=t\, 
varies with increasing values of nearest neighbor pair-hopping amplitude Y\ and has a maximum at Y\ w 1.33/Ji, instead 
of having a constant value. Second, the value of this maximum U max « — 0.33z]ii is independent of X\. Furthermore, 
the stability regions for all values-of X\ jt\ derived with the present method are always larger than the corresponding 
ones predicted by Montorsi et alx3 and de Boer et alS3 and exist for all positive values of Y\. The relation Y\ — —2V\ 
is also satisfied because the ^-state has to be an exact eigenstate of the model. Therefore, the positivity of Y\ implies 
that the nearest neighbor density-density type interaction V\ has to be always negative, i.e. attractive, in order to 
find an 77^ ground state. 

As the 77-pairing states consist of local pairs of electrons, it is of interest to investigate the effect of on-site Coulomb 
repulsion, characterized by U, on the stability of these pairs. To do this the V 2 -Vi plane is chosen at various values 
of U. As can be seen in Fig. I the local pairs are stable even in the presence of relatively large positive values of U. 
This requires, however, an attraction in the nearest neighbor density-density interaction channel. It should also be 
noted that the same type of interaction between next-to-nearest neighbors, denoted by V2, can either be attractive or 
moderately repulsive. These findings lead us to the conclusion that the 77-pairing state of momentum P = ir remains 
the ground states of Eq. (0) even for positive values of the on-site Coulomb interaction. Hence superconductivity can 
exist in the extended Hubbard model with local repulsion (U >0), if a sufficiently strong nearest neighbor attraction 
(Vi < 0) is present. 

The inclusion of next-to-nearest neighbor interactions increases remarkably the number of model parameters and 
hence the number of possible cuts of the parameter space. Therefore, we illustrate only some special, overall effects of 
these interactions. In order to model real systems all next-to-nearest neighbor interactions are chosen to be smaller in 
magnitude than the corresponding nearest neighbor ones. Nevertheless, the ratio of nearest to next-to-nearest neighbor 
interactions can be very different - it depends on the material. In Fig. ^ three plots are shown for different values 
of the couplings. These belong to the case of D = 2. The corresponding 3-dimensional plots display qualitatively 
the same features, except for the parameters of Fig. ^c. The quantitative discrepancy between the plots taken in 
D = 2 and in D = 3 is due to the fact that the number of next-to- nearest neighbors z 2 is much larger in D = 3 spatial 
dimensions than in D = 2 on a hypercubic lattice. This suggests that in the framework of the present model the effects 
of next-to-nearest neighbor interactions are stronger in higher dimensions. 

In Fig. |^a the stability region of the 77-pairing state of momentum P — is shown for a certain set of model 
parameters in the absence (solid line) and in the presence (dotted line) of next-to-nearest neighbor couplings. One 
can notice the expansion of the stability region due to the presence of next-to-nearest neighbor interactions. Since a 
large value of |ii| (in the presence of a fixed value of the nearest neighbor pair hopping Y\) favors the hopping of single 
electrons over the hopping of electron pairs, large values of \t\\ give rise to the breaking of local pairs. This means 
that the number of doubly occupied sites is not conserved any longer and, as a result, the 770-state ceases to be the 
ground state of -ffgiob- In Fig. J^b the effects of next-to-nearest neighbor interactions on the stability of the T/V-state 
are shown. In contrast to the situation depicted in Fig. ^|a no shrinking of the stability domain with increasing \t\\ 
can be observed. This can be explained with the different internal structure of the 77^ pairs. In Fig. ||c we show a 
situation, where next-to-nearest neighbor couplings can either extend or shrink the stability region of the 77^-state. 
As mentioned earlier, the dimension of the lattice plays a crucial role here. In D = 2 a huge portion of the U-Y\ plane 
phase diagram is occupied by the T^-state for any ratio of X\/t\. In D = 3, however, it was found for a wide parameter 
region that the rj^-state is the ground state of the model Hamiltonian i? g i b only for the special case of X\/t\ = \. 

B. The fully polarized ferromagnetic state 

Let us consider now the fully (z-)polarized ferromagnetic state as the trial wave function defined as 
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L 



0- 



Calculating the commutator of F with i? g i b one can see that this state is an exact eigenstate of the global Hamiltonian 
for any values of the interaction parameters. The trial wavefunction |\&fm) yields in the thermodynamic limit at half- 
filling the upper bound 



upper 



1 1 2 

--UL + -LJ2^4 1) - l*L (21) 



1=1 



for the ground state energy. 

In what follows, we reveal under what conditions Eq. ( p0| ) is the ground state of Hgiob- For the sake of simplicity 
we concentrate on a fixed set of numerical values of nearest neighbor couplings, as it has already been estimated by 
HubbardES for electrons in rf-bands of transition metals. The values of next-to- nearest neighbor interactions are chosen 
to be fractions of the corresponding nearest neighbor ones. The ratio of nearest to next-to- nearest neighbor couplings 
is set to be about 5-8. We think this rangej-ef ratio to be appropriate for a wide class of materials. Furthermore, this is 
in agreement with the work of Appel et alEB who also made quantitative predictions regarding the strength of nearest 
and next-to- nearest neighbor correlated hopping X\ and X2, respectively. Further calculations at different sets of 
model parameters have also shown that the phase diagrams plotted in Figs. [?], || and || are generic. This suggests that 
the above choice of model parameters captures the essential physics. 

In Fig. [7] we present the changes inpthe stability domain induced by using plaquette Hamiltonians. For comparison 
with the corresponding result of Rcf£3 with bond Hamiltonians, the plaquette Hamiltonians contained no next-to- 
nearest neighbor interactions. The shaded region shows the enlargement in the stability domain of the fully saturated 
ferromagnetic state. As can be seen from the figure, there is a reasonable extension with respect to t\. While in 
the case of bond Hamiltonians the value of the correlated hopping X\ should be very close in magnitude to ii in 
order to reach the boundary of the stability region at U m i n w 4 eV , we have a much broader region for that using 
plaquette Hamiltonians. The broadening implies that the additionally incorporated spatial correlations really lead to 
the stabilization of the ordered phase, in our case the ferromagnetism. 

The global Hamiltonian -ffgiob containing purely nearest neighbor interactions can be transformed into an effective 
Heisenberg model in the large- U limit at half filling (see e.g. Kollar et alxB and references therein) with the effective 
exchange coupling of 

J efi = I ( 1 - £iY + J, . (22) 



U V h 

This favors ferromagnetism for all J c g < 0. Neglecting in Eq. (|t]) all intersite interactions but the nearest neighbor 
exchange interaction, the latter equation suggests a simple perturbative criterion for the stability of ferromagnetism 
in the large- U limit; ferromagnetism is favored over antiferromagnetism for all U > U c — t\/\Ji\ (note that in Eq. 
(0) Ji = — \ <0 means the ferromagnetic coupling). It is known that the OGS approach using bond Hamiltonians 
gives the criterion U/z\ > U c as the stability requirement in the same regime. Therefore it possibly underestimates 
the stability of the fully polarized ferromagnetic state. Taking into account larger local Hamiltonians defined on 
elementary plaquettes of the hypercubic lattice one gets the criterion U/zi > U c /2 for the stability, which means 
more extended stability domain in the large- U limit. Furthermore, it suggests that (at least in the large- U limit) 
the stability criterion has the form of U > ^-U c , where b is the function of the size of the cluster on which the local 
Hamiltonian is defined,— The possible scaling behavior of the stability criterion and the concrete form of ^(A^iuster) 
are discussed elsewhere^. 

Since the fully polarized ferromagnetic state at half filling is an exact eigenstate of H g \ b we have got no a priori 
restrictions for the values of the interaction parameters. The extensive calculations, however, lead to a simple restric- 
tion between jiy and ■ In order to have a ferromagnetic ground state of a model containing spin interactions 
which vary continuously from a Heisenberg- type interaction to a simple Ising-type one, the — 1 < A^ Y < 1 requirement 
must hold. This means that the restriction 

-\JP\<-J$<\JP\ (23) 

has to be always satisfied. In Fig. [8] the consequence of Eq. ( p3| ) is illustrated in the U-t\ cut of the parameter space 

at Af = 1 and jj = — | Ji | < for various values of Af- Y . The size of the stability region is maximal at A* Y = 1 
and gradually decreases as A* Y reaches A* Y = — 1. Any further decrease of A* Y yields that our ferromagnetic state 
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which is fully z-polarized is no longer the ground state of -f/giob! for Jxy = jj 1 ^ -ffgiob is SU(2) symmetric which implies 
the degeneracy of the ground state with respect to SU(2) rotations. For anisotropic exchange couplings favouring the 
xy-plane, the ground state polarization may still be macroscopic, however, it no longer points in the z-direction. This 
means that the fully z-polarized ferromagnetic trial wavefunction becomes unstable. We note that Eq. (|2^) must hold 
even in the presence of next-to-nearest neighbor couplings. 

In Fig. |]the effects of nonzero next-to- nearest neighbor hopping £2 are illustrated for a fixed set of model parameters, 
first in the U-t\ plane [plot (a), in units of eV\ and second in the U-t% plane [plot (b), in units of t%\. In Fig. [|a the 
dotted line represents the phase boundary in the absence of next-to-nearest neighbor interactions while solid, long- 
dashed, dashed and dotted-dashed lines correspond to phase boundaries in the presence of next-to-nearest neighbor 
interactions for various values of £2. As can be seen from the plot, next-to- nearest neighbor interactions can help in 
stabilizing ferromagnetism for the chosen set of model parameters as long as tijt\ has a small, positive value. For 
negative or large positive values of £2/^1 > however, the stability domain reduces significantly and stronger Coulomb 
repulsion is needed for the stabilization. This means e.g. that for a reasonably narrow band (ti «0.4 eV) with the 
ratio of t 2 /^i =0.1, the required minimal stabilizing Coulomb repulsion is about U m - m w 3 eV. The same value of 
U at ^Ai = —0.25 is about [/,„;„ « 30 eV which is a magnitude larger. Nevertheless, the above range of Coulomb 
interactions can be considered reasonable for real materials. 

In Fig. |^b the effects of change in non-interacting dispersion due to the inclusion of next-to-nearest neighbor hopping 
are shown for a square (solid line) and for a cubic (dotted line) lattice. The stability domains of ferromagnetism depend 
on the dimensionality of the lattice and do not coincide. It is interesting to note that the most favorable values of £2/^1 
for which the Coulomb interaction takes its minimal value /7 m ; n also depend on the spatial dimension. The shapes of 
the stability domains are also of interest: in a well-defined region of £2/^1 U m i n changes only by a slight amount but 
as soon as the edges of this region are reached U increases drastically. This feature sugges±s_that inside the stability 
regions a nonzero next-to-nearest neighbor hopping via the asymmetric density of statest3E3 helps (in the presence 
of other next-to-nearest neighbor interactions) in stabilizing ferromagnetism but outside it destabilizes ferromagnetic 
ordering. The edges are determined mainly by the dispersions (hence the shape of the particular density of states) 
and tuned further by other interactions being present in H g i b- 

It is also known that the inclusion of nearest^ighbor ferromagnetic exchange interaction in the pure Hubbard model 
favors the parallel ordering of electron spinsoM'tJ. In Fig. |l^ we considered the pure Hubbard model supplemented 
with next-to-nearest neighbor hopping £2 and nearest and next-to-nearest neighbor exchange interactions Ji and J2, 
respectively. All the other type of couplings are turned off. As can be seen in the figure, the stability domain of 
ferromagnetism extends with increasing value of the Coulomb repulsion and in the limit of U — > 00 fills the whole 
J\ , J2 < quarter of the phase diagram. The fully polarized ferromagnetic state remains the ground state of Eq. (Q) 
for finite values of U only in the presence of finite values of J\ and J2. Figure [l^ also shows that the required minimal 
values of | J2 1 are about an order of magnitude less than the required minimal values of \J\\, and J2 should also be 
ferromagnetic in nature, i.e. J2 <0. 

V. CONCLUSIONS 

In the present paper we have studied in the thermodynamic limit the ground state phase diagram of the Hubbard 
model supplemented by nearest and next-to-nearest neighbor interactions. The purpose of the study was to clarify 
to what extent and in which way the inclusion of additional spatial correlations changes the stability of physically 
interesting states, the 77-pairing state of momentum P — and P = n, or the fully z-polarized ferromagnetic state. 
The phase boundaries are extracted from the equality of an upper and a lower bound of the ground state energy, 
hence these are exact. The additional spatial correlations are introduced via the computation of the lower bound on 
elementary plaquettes, instead of bonds, of the D dimensional hypercubic lattice. Except the case of 77-pairing state 
of momentum P = the exact phase boundaries cannot be given in closed, analytic forms. Instead, they are shown 
graphically in special cuts of the parameter space of the model under study. 

The phase boundaries presented are sufficient phase boundaries. This means that outside the region defined by the 
exact conditions a certain state might remain the ground state of the model. Diagonalizing local Hamiltonians defined 
on larger clusters of the lattice the phase boundaries might be further improved. This improvement with increasing 
cluster size addresses a further issue: does the stability domain of a given ordered phase extend further with taking 
larger and larger clusters, or is there a convergence regarding the location of its phase boundary? If the latter holds, 
we could determine the phase boundaries with a simple extrapolation even in the limit of H = Cluster • Based on 
preliminary results we believe that the further extension of stability domains decreases rapidly with increasing cluster 
size. For instance, computing the lower bounds of the ground state energy with diagonalizing local Hamiltonians 
defined on clusters of 6 lattice sites, the further expansion of the stability regions is only about a few percentage, 







generally 4-5 % or less. 

Considering the effects of more spatial correlations, which was equivalent in our case with th a ^hoic e of plaquette 
Hamiltonians, we have improved significantly the previously derived rigorous results of RefsMBEiH This means 
e.g. for the fully polarized ferromagnetic state that the minimal value of the Coulomb repulsion required to stabilize 
ferromagnetism along reasonable values of nearest neighbor interactions is predicted to be about 6-10 eV (in D = 3) 
in a relatively broad range of nearest neighbor hopping (see Fig. |?]). 

Another goal of the present study was to determine the overall effects of next-to-nearest neighbor interactions on 
the stability domains. The inclusion of next-to- nearest neighbor interactions can be done naturally using plaquette 
Hamiltonians. Up to our knowledge, the effects of next-to-nearest neighbor interactions, except those of next-to-nearest 
neighbor hopping, have not yet been considered rigorously in the literature. The relative strengths of next-to-nearest 
neighbor interactions are much smaller than that of nearest neighbor interactions, the stability conditions, however, 
strongly depend on them. Their effect in various sectors of the phase diagram is different and they result either an 
extension or a shrinking of the stability domain. For instance, taking the 77-pairing state of momentum P = tt the 
possible maximal value of the Coulomb repulsion C/ max , up to which the ^-state remains the ground state of the 
extended Hubbard model (i.e. the model has a superconducting ground state), is increased from 6 eV to 8-9 eV (in 
£> = 3) by the inclusion of relatively small next-to- nearest neighbor interactions (Fig. |^b). 

It is also known that next-to-nearest neighbor hopping of single particles, which is characterized by the hopping 
amplitude t2, is of importance in real materials. Our results are in good agreement with this fact. We showed that ti 
has a characteristic effect e.g. on the stability of the fully saturated ferromagnetic state. Even a small ratio of t 2 /t 1 , 
i.e. a small amount of frustration in the dispersion, introduces a qualitative change into the phase diagram. The 
change is mostly a shrinking of the stability domain, however, for small ratios it^/tx <0.15) the presence of ti helps 
in stabilizing ferijG«nagnetism (Fig. ^). This is in good agreement with recent DMRG studies taken on 1-dimensional 
triangular latticeE3. It is interesting to note that in our calculations the extension of ferromagnetic domain occurs 
always at positive ratios of tilt\ for fixed values of the other parameters of the model. 

The Hubbard model supplemented only by exchange interactions J\ and J2 has also been investigated. Our results 
are in good agreement with RefEfl, i.e. the critical values of nearest and next-to-nearest exchange interactions to give 
rise to ferromagnetism approach zero as U — >oo in the case of a half-filled band in any dimensions. However, at finite 
values of the Coloumb repulsion J\ and J2 should also be finite, if the ground state is the fully polarized ferromagnetic 
state. 

In summary, we have established a simple method which allows us to incorporate and to treat the effects of next-to- 
nearest neighbor correlations and interactions in an exact fashion. We showed that the ground state of the extended 
Hubbard model in the thermodynammic limit at half filling is superconducting or ferromagnetic, depending on the 
interaction strengths. The improved phase boundaries for certain sets of model parameters have also been constructed. 
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FIG. 1. An elementary plaquette of the D = 2 dimensional square lattice. Solid lines represent different types of nearest 
neighbor couplings while dashed lines symbolize the next-to-nearest neighbor ones. The local Hamiltonian is defined on this 
cluster. 
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FIG. 2. Exact stability region of ^-pairing state of momentum P — (r/o) at half filling. The shaded region represents the 
enlargement of the stability domain due to the choice of local Hamiltonian defined on elementary plaquettes (see Fig. |l]) of the 
lattice. All the next-to-nearest neighbor interactions are kept zero. 
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FIG. 3. Exact stability domains of ^-pairing state of momentum P — tt (rj w ) at half filling for two different sets of nearest 
neighbor couplings in the absence of next-to-nearest neighbor interactions [plot (a): X\ =0.1 eV, Vi = —2 eV and Ji =0.05 eV; 



plot (b): Xi=0.1 eV, V^- 



-2 eV, J$ 



-0.02 eV and J, 



(i) . 



-0.015 eV]. The shaded regions represent the enlargement of 



the stability domain due to the choice of local Hamiltonian defined on plaquettes instead of bonds. 
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FIG. 4. Effects of correlated hopping Xi and pair hopping Y\ on the stability of 77-pairing state of momentum P — n (r/^). 
Each next-to-nearest neighbor interaction is turned off. We note, that beyond a well-defined value of Yi the phase boundaries 
for the cases X\ =fct\ and Xi =ti coincide. 
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FIG. 5. Stability domains of ^-pairing state of momentum P — ir (r/^) in D — 3 in the presence of various on-site Coulomb 
repulsions and intersite density-density type interactions (Vi,V2). The remaining parameters of Eq. (Q) are chosen as follows: 
<2 = — Jti and X i = Ji = J2 = 0. Note, that X2 ,Yi and I2 are uniquely fixed via the rigorous restrictions derived earlier in order 
that Eq. (113) be an exact eigenstate. In D — 2 similar stability regions emerge. 
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FIG. 6. Effects of next-to-nearest neighbor (NNN) couplings on the stability of ^-pairing states of momentum P on a 
two dimensional square lattice. The values of the interaction constants are fixed as follows; plot (a): P = 0, Vi = —3 eV, 
J xy = io eV > J - (1> = 55 eV ' *a = V 2 = \Vi and A 2) = ijP (a=xy,z); plot (b): P = it, Vi = -2 eV, Xi = 0.5 eV, 
jiy>=-± eV, Jz (1) = -^ eV, *2 = -§ti, Va = |Vi and ji 2) = iji 15 (a = xy,z); plot (c): P = tt, Xi = |ti, Ji = -2Y U t 2 = -$ti, 
Y 2 = \Y 1 and J 2 = |Ji. 
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FIG. 7. Exact stability region for the fully saturated ferromagnetic state (FM) at half filling onafl dimensional hypercubic 
lattice for a certain set of model parameters X\ — 0.5 eV, Vi = 2 eV and Y\ = — Ji = ^ eV in the absence of next-to-nearest 
neighbor interactions. The shaded region represents the extension of the stability of the fully saturated ferromagnetic state as 
the ground state of Eq. ((?]) due to the choice of local Hamiltonian defined on elementary plaquettes, instead of bonds, of the 
lattice. 
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FIG. 8. Stability regions of the fully saturated ferromagnetic state (FM) for various values of A xy 
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FIG. 9. Effects of next-to-nearest neighbor couplings on the stability of the fully saturated ferromagnetic state (FM) at half 
filling in the t%- U and £2- U planes, plot (a) and (b), respectively. Plot (a) shows the phase boundaries for various ratios of tz/ti 
choosing the numerical values of nearest neighbor interactions as in Fig. Next-to-nearest neighbor interactions for the plot 
are as follows: X2 =0.08 eV, Vz = gVi, Yi = \Yi and J2 = |Ji. In plot (b) all the interaction constants are expressed in units 
of ti instead of units of e V. Above each line the ground state is the fully polarized ferromagnetic state. 
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FIG. 10. Stability of the fully saturated ferromagnetic state (FM) in the presence of nearest (Ji) and next-to-nearest (J2) 
neighbor exchange coupling at ti = — j^ti. All the other interactions are turned off, i.e. Xi =Xi = Vi = Vi = Yi = Yi = 0. 
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